function [rmin,rmax,etaDraw] = approximateBounds(restr,phi,opt)
% Approximate bounds of the conditional identified set via simulation.
% Inputs:
% - restr: structure containing information about restrictions
% - phi: structure containing reduced-form VAR parameters
% - opt: structure containing model information and options

H = opt.H;
ivar = opt.ivar;
jshock = opt.jshock;
qDraws = opt.qDraws;
vma = phi.vma;
cumIR = opt.cumIR;

% Draw from space of orthonormal matrices satisfying the restrictions 
% qDraws times.
Q0 = drawQs_ar(restr,phi,opt);
% Keep columns of Q needed to compute impulse responses.
Q0 = Q0(:,jshock,:); 

% Compute impulse responses for each draw of Q.
etaDraw = zeros(H+1,length(ivar),length(jshock),qDraws);

for hh = 1:H+1 % For each horizon

    % Extract required rows of horizon-h VMA coefficient matrix.
    Cphi = vma(ivar,:,hh);
    % Multiply by relevant columns of Q0 (for each draw of Q0) to obtain
    % impulse response.
    etaDraw(hh,:,:,:) = pagemtimes(Cphi,Q0);

end
 
if ~isempty(cumIR)
    
    % Compute cumulative impulse responses for relevant variables.
    etaDraw(:,cumIR,:,:) = cumsum(etaDraw(:,cumIR,:,:),1);

end

% Compute minimum and maximum impulse response over draws of Q.
rmin = min(etaDraw,[],4);
rmax = max(etaDraw,[],4);

end

